function [Effect,effect_var,invtausq_effect,llik_trace] = var_logit(Y,X,BURNIN,NEPOCH,Effect,invtausq_effect)
%
% [Effect,EffectVar,invTauSq,llik] = update_logit_model(Y,X,BURNIN,NEPOCH,EFFECT,invTauSq)
%
% Inference of prior logit model
%
% Y(i,:) ~ logit(X(i,:)*Effect')
%
% Y           = sample x phenotype probabilistic/binary matrix
% X           = sample x feature data matrix (real number)
% BURNIN      = burn-in optimization time without TauSq
% NEPOCH      = optimization time with TauSq
% EFFECT      = pheno x feature effect size
% EFFECTVAR   = pheno x feature effect size variance
% invTauSq    = feature x 1 shrinkage (feature selection)
%
% code: Yongjin Park, ypp@csail.mit.edu
%

    TOL = 1e-4;

    % ================================================================

    [Nsample,Npheno] = size(Y);
    [nsample,Nfeat] = size(X);
    assert(Nsample == nsample);

    % initialize
    if nargin > 5,
        assert(Npheno == size(Effect,1) && Nfeat == size(Effect,2));
        assert(size(invtausq_effect,1) == Nfeat);
    else
        Effect = zeros(Npheno,Nfeat,'single');
        invtausq_effect = 0.01*Npheno*ones(Nfeat,1,'single');
    end
    
    llik_trace = NaN(BURNIN+NEPOCH,1,'single');

    for iter = 1:(BURNIN+NEPOCH),

        % construct local quadratic approximation
        F = X * Effect';
        P = 1./(1+exp(-F));
        P(P < 1e-4) = 1e-4;
        P(P > 1 - 1e-4) = 1 - 1e-4;

        W = P.*(1-P);
        workRes = F + (Y - P) ./ W;

        % coordinate-descent step
        for jj = 1:Nfeat,

            wh_j = bsxfun(@times, W, X(:,jj));
            wh2_j = bsxfun(@times, W, X(:,jj).^2);

            denom_effect = sum(wh2_j)' + invtausq_effect(jj);
            num_effect = sum(wh_j .* (workRes - F + X(:,jj)*Effect(:,jj)'))';
            aa = num_effect ./ denom_effect;

            % update F
            F = F + X(:,jj)*(aa - Effect(:,jj))';

            Effect(:,jj) = aa;
        end
        
        llik = sum(sum(Y .* log(P) + (1-Y) .* log(1-P)));
        llik = llik - 0.5 * sum(sum((Effect.^2)*invtausq_effect));
        llik_trace(iter) = llik;

        if iter > BURNIN,
            invtausq_effect = Npheno./sum(Effect.^2,1)';
            invtausq_effect = arrayfun(@(x) max(0.01/Npheno,min(100*Npheno,x)), invtausq_effect);
        end

        if iter > BURNIN + 5,
            llik_prev = mean(llik_trace((iter-5):(iter-1)));
            llik_prev_std = 1e-10 + std(llik_trace((iter-5):(iter-1)));
            if abs(llik_prev - llik)/llik_prev_std < TOL,
                llik_trace = llik_trace(1:iter);
                fprintf(2,'Converged\r');
                effect_var = 1./(W'*(X.^2) + repmat(invtausq_effect',Npheno,1));
                break;
            end
        end
        
        fprintf(2,'Iter = %03d, LLIK = %.4e\r',iter,llik);
    end

end
